Experimental OPF Implementation using MATPOWER Input by lukelowry · Pull Request #594 · ORNL/GridKit

lukelowry · 2026-09-28T17:25:11Z

Description

Experimental @superwhiskers @nkoukpaizan which uses MATPOWER input which makes interfacing with PhasorDynamics harder but most OPF data comes in MATPOWER, so this is probably most practical choice

Application OptimalDispatch

Checklist

All tests pass.
Code compiles cleanly with flags -Wall -Wpedantic -Wconversion -Wextra.
The new code follows GridKit™ style guidelines.
There are unit tests for the new code.
The new code is documented.
The feature branch is rebased with respect to the target branch.
The CHANGELOG.md has been updated to reflect the changes. If this is a minor PR that is part of a larger fix already included in the file, state so.

Further comments

Needs some cleaning first

pelesh · 2026-09-28T18:22:04Z

A few comments:

If you have JSON input parser already available, I suggest you use that instead of MATPOWER format. I was thinking you could use some of the existing power flow implementation, but since you took a different route, I suggest you ditch MATPOWER and go straight to JSON.
I would have optimal power flow app in EconomicDispatch subdirectory in examples directory, at the same level as PhasorDynamics.
I would rename Optimization/OptimizationProblem.*pp files (and the class within) into Optimization/SteadyState.*pp. We will likely have Multiperiod and MixedInteger as new optimization modules in the future.

lukelowry · 2026-09-28T20:04:42Z

@pelesh I changed it so it uses the PowerFlow parser with very few changes. I do believe its valuable to have that as the input format only because it makes validation easier and there is practically unlimited data in .m, and OPF data does not always come with dynamic case files.

I believed I addressed your other issues too. I am able to use the example app to change the objective/model, solve for new operating point, and start simulation in steady state at that new point!

pelesh

This is very nice work. I suggest limiting the scope of this PR to providing ACOPF module only. Imho, the integration with dynamic simulations modules should be done after ACOPF is tested and input file format specified.

superwhiskers

tests pass on my end. i think it works quite well, but there are a few implementation questions we should sort out before this moves forward.

additionally, if we are to continue with matpower, are we intending to handle capability curves for generators or the convention of a dispatchable load as a generator with negative real-power injection? if they are outside the intended scope, could we reject them and document that limitation?

i think these conventions, particularly the latter, strengthen the argument for going with our own format, as using an existing format such as matpower requires us to adhere to its conventions, which may not be desirable.

superwhiskers · 2026-10-16T15:32:11Z

+/// Relative tolerance of the optimal cost against a reference cost
+constexpr double COST_TOLERANCE = 1.0e-8;


is there not a situation in which one may want to configure this?

Yes, the user will want to specify the tolerance. And this tolerance is too tight.

superwhiskers · 2026-10-16T15:39:33Z

+      /**
+       * @brief Coordinate entries in global indices
+       *
+       * Filled by `mapped_store` and `lower_store`. Assembly sums duplicate
+       * coordinates.
+       */
+      struct CooEntries
+      {
+        std::vector<size_t> rows;
+        std::vector<size_t> cols;
+        std::vector<double> values;
+
+        void clear()
+        {
+          rows.clear();
+          cols.clear();
+          values.clear();
+        }


you should parameterize this by index type / scalar type so that you don't need the static assertion that IdxT/ScalarT are size_t/double in the evaluators

superwhiskers · 2026-10-16T15:40:19Z

+        }
+
+      private:
+        __attribute__((always_inline)) inline static ScalarT lagrangian(const ModelT*  model,


prefer [[gnu::always_inline]]

superwhiskers · 2026-10-16T15:43:46Z

+       * @tparam ModelT - model type with `VARIABLE_SIZE` local variables and
+       * `CONSTRAINT_SIZE` local constraints


VARIABLE_SIZE may be more accurately named VARIABLE_COUNT, NUM_VARIABLES, VARIABLE_ARRAY_SIZE, or something else. VARIABLE_SIZE to me read as "bit/byte width of individual variables"

superwhiskers · 2026-10-16T15:47:21Z

+       * @tparam ModelT - model type with `VARIABLE_SIZE` local variables and
+       * `CONSTRAINT_SIZE` local constraints


following what i did here, i'd recommend establishing these requirements on the ModelT parameter as a concept to reduce the amount of error noise for the end user

superwhiskers · 2026-10-16T15:54:25Z

+       * @param[in] col_indices - Global column of each local column
+       * @param[in,out] entries - Entries in global indices
+       */
+      [[maybe_unused]] __attribute__((always_inline)) static void mapped_store(double        val,


Suggested change

[[maybe_unused]] __attribute__((always_inline)) static void mapped_store(double val,

[[maybe_unused, gnu::always_inline]] static void mapped_store(double val,

superwhiskers · 2026-10-16T16:00:59Z

+      static constexpr IdxT INTERNAL_VARIABLE_SIZE   = 0;
+      /// Internal constraint functions: squared apparent power at bus 1 and bus 2
+      static constexpr IdxT INTERNAL_CONSTRAINT_SIZE = 2;
+      static constexpr bool CONSTANT                 = false;


i would like you to document the meaning of this parameter somewhere. one place to do that would be on the concept that i proposed in another review comment for compile-time validation of optimal power flow components

superwhiskers · 2026-10-16T16:03:17Z

+    /**
+     * @brief Line or off-nominal transformer between two buses
+     *
+     * Uses the admittance of `PhasorDynamics::Branch`. Power is positive


"uses the admittance of PhasorDynamics::Branch" and similar comments referencing how the modeling code is based off that in the phasor dynamics modeling library are not worth retaining as they could diverge in the future, imo

superwhiskers · 2026-10-16T19:05:55Z

+          const auto entry = actual.find(coordinate);
+          if (entry == actual.end())
+          {
+            std::cerr << "Entry (" << coordinate.first << ", " << coordinate.second << ") is missing\n";


you should be using the logger here and elsewhere in these tests

superwhiskers · 2026-10-16T19:44:25Z

+
+          auto& generator                                 = data.generator.emplace_back();
+          generator.id                                    = std::format("gen_{}", i + 1);
+          generator.in_service                            = row.status > 0;


because SystemModelData parameterizes the power flow data types with size_t, -1 (which is a valid status value) will be converted to a positive value, which will misinterpret the model.

possible solutions:

don't use matpower, use our own json modeling capabilities which could represent it as a boolean

change how status is represented, don't tie it to IdxT, probably parse it to a boolean but read it as int or something

nkoukpaizan

Overall nice work! The component implementations are clear and I like the consolidated ComponentModelImpl.hpp.

A few comments:

For Enzyme, can we re-use some of the existing boilerplate for LowerSparseStorage and use buffers instead of CooEntries wrapping std::vector)? Feel free to amend the wrappers.
Can we bind GridKit::LinearAlgebra::Vector at the component level instead of std::vector? This could be a follow-up PR.
Several methods, e.g., in OptimizationEvaluator* and SteadyState* don't have a Doxygen documentation.

nkoukpaizan · 2026-10-16T22:16:25Z

+      static constexpr RealT POWER_TOL    = 2 * std::numeric_limits<RealT>::epsilon();
+      static constexpr RealT JACOBIAN_TOL = 2 * std::numeric_limits<RealT>::epsilon();
+      static constexpr RealT HESSIAN_TOL  = 9 * std::numeric_limits<RealT>::epsilon();


I would round to 10 * and use the same tolerance.

nkoukpaizan · 2026-10-16T22:18:26Z

+if(GRIDKIT_ENABLE_ENZYME)
+  add_subdirectory(OptimalPowerFlow)
+endif()


I would move the guard into the OptimalPowerFlow directory.

nkoukpaizan · 2026-10-16T22:22:21Z

+      static constexpr RealT DERIVATIVE_TOL = 8 * std::numeric_limits<RealT>::epsilon();
+      static constexpr RealT BALANCE_TOL    = 1.0e-12;


Same here, I would round 8 to 10. And replace 1e-12 with a multiple of std::numeric_limits<RealT>::epsilon();

nkoukpaizan · 2026-10-16T22:25:01Z

+      bool isZero(double value)
+      {
+        return !(value < ZERO<double>) && !(value > ZERO<double>);
+      }


We've been trying to avoid these. If something needs == 0/!= 0, please highlight and we can evaluate.

nkoukpaizan · 2026-10-16T22:27:41Z

+          return true;
+        }
+        std::cerr << std::setprecision(std::numeric_limits<RealT>::max_digits10) << value << " differs from "
+                  << reference << " by " << std::abs(value - reference) / (ONE<RealT> + std::abs(reference)) / std::numeric_limits<RealT>::epsilon()


I wouldn't re-implement Testing::isEqual.

nkoukpaizan · 2026-10-16T22:30:48Z

+        const RealT* x = system.x().getData();
+
+        success *= isEqual(x[1], ZERO<RealT>, BALANCE_TOL);
+        success *= std::abs(x[3]) > BALANCE_TOL && std::abs(x[5]) > BALANCE_TOL;


Is checking the absolute values really needed instead of isEqual?

nkoukpaizan · 2026-10-16T22:48:43Z

    is >> br.fbus >> br.tbus >> br.r >> br.x >> br.b >> br.rateA >> br.rateB
        >> br.rateC >> br.ratio >> br.angle >> br.status >> br.angmin
        >> br.angmax;
-    checkEndOfMatrixRow(is);


Removing checkEndOfMatrixRow partially addresses #587, but we still need to ensure the read was successful (e.g., row not too short).

pelesh

This looks very good and has clean component model structure. Main issues that I would raise are:

Separating equality and inequality constraint evaluation methods.
Selection of coordinates polar vs. Cartesian.
Inheritance patterns should mirror PhasorDynamics module. If we want to change it, then the change should made be across the GridKit modules.

pelesh · 2026-10-16T23:55:04Z

+       * Filled by `mapped_store` and `lower_store`. Assembly sums duplicate
+       * coordinates.
+       */
+      struct CooEntries


What is the purpose of this struct? This needs to be documented.

pelesh · 2026-10-17T00:05:05Z

+      using ScalarT    = scalar_type;
+      using IdxT       = index_type;
+      using RealT      = typename ComponentModelT::RealT;
+      using ModelDataT = BranchData<RealT, IdxT>;


These should not be public. ScalarT, IdxT, etc. are local type aliases.

pelesh · 2026-10-17T00:07:51Z

+      g[0] = p1 * p1 + q1 * q1;
+      g[1] = p2 * p2 + q2 * q2;
+      g[2] = p1;
+      g[3] = q1;
+      g[4] = p2;
+      g[5] = q2;


Need to specify here what are equality and what are inequality constraints.

pelesh · 2026-10-17T00:17:31Z

+     * @tparam model_type - Derived component model
+     */
+    template <typename model_type, typename scalar_type, typename index_type>
+    class ComponentModel : public Component<scalar_type, index_type>
+    {
+    public:


I suggest removing Component class and using the same inheritance structure as in PhasorDynamics module.

pelesh · 2026-10-17T00:23:34Z

+      const auto& terminals = component.terminals();
+      for (IdxT t = 0; t < terminals.size(); ++t)
+      {


I believe that we decided to use term "ports" rather than "terminals". It would be good to use consistent nomenclature across all GridKit modules.

pelesh · 2026-10-17T00:32:11Z

+                                                               Number* x_u,
+                                                               Index   m,


These are deprecated Ipopt types. Use ipnumber and ipindex instead.

pelesh · 2026-10-17T00:36:34Z

+    template <typename scalar_type, typename index_type>
+    class SteadyState : public Ipopt::TNLP
+    {


This is not SteadyState but rather Ipopt specific interface. Perhaps use SteadyState for the namespace and IpoptInterface for the class name.

pelesh · 2026-10-17T00:48:55Z

+    public:
+      using ScalarT     = scalar_type;
+      using IdxT        = index_type;
+      using RealT       = typename ScalarTraits<ScalarT>::RealT;
+      using CooEntriesT = Enzyme::Sparse::CooEntries;


These type aliases should not be public.

pelesh · 2026-10-17T00:53:00Z

+      std::vector<std::unique_ptr<ComponentT>> components_;
+
+      /// Buses by number
+      std::map<IdxT, BusT*> buses_;
+
+      /// Map from component entries to CSR entries
+      std::vector<IdxT> jacobian_map_;
+      std::vector<IdxT> hessian_map_;


I would avoid using STL containers for system level objects. If it is a temporary solution, it should be marked as such.

pelesh · 2026-10-17T01:04:21Z

+      virtual int evaluateConstraints() = 0;
+      virtual int evaluateJacobian()    = 0;


It would be good to separate equality and inequality constraints, as well as their respective Jacobians.

lukelowry · 2026-10-18T13:40:07Z

Thank you all for detailed review, I believe I have addressed all except a distinct few. I will go response by response soon, but the overall design is overall more consistent with PhasorDynamics except for the inheritance structure. I am considering removing the larger ACTIVSg2000 case until after this is merged, please let me know what you think.

lukelowry requested review from nkoukpaizan and superwhiskers September 28, 2026 17:25

pelesh self-requested a review September 28, 2026 18:13

pelesh requested a review from abirchfield September 28, 2026 18:22

lukelowry force-pushed the lukel/OPF-NEW-dev branch 4 times, most recently from 95913b8 to 5c4fb92 Compare October 10, 2026 21:02

lukelowry force-pushed the lukel/OPF-NEW-dev branch from 5c4fb92 to 49ad14b Compare October 13, 2026 22:04

pelesh reviewed Oct 14, 2026

View reviewed changes

lukelowry added 9 commits October 14, 2026 23:54

Initial OPF implementation

860af60

Minimize Churn on PhasorDynamics

81ba74e

clenaup matpower for both models

dd18dcd

rename to EconomicDispatch and downgrade to an example app

51996ad

ipopt cmake settings

2b40e6b

remove integration with PhasorDynamics Apps

583d1af

cleanup uncessary tests

d428593

BUG: Application symlink issue

b74ff99

Enhance and cleanup

a6e3ae1

lukelowry force-pushed the lukel/OPF-NEW-dev branch from 979e671 to a6e3ae1 Compare October 15, 2026 07:09

platofrm dependent correction

d5c7d68

superwhiskers requested changes Oct 16, 2026

View reviewed changes

nkoukpaizan reviewed Oct 16, 2026

View reviewed changes

nkoukpaizan assigned lukelowry Oct 16, 2026

nkoukpaizan added enhancement New feature or request new model new solver labels Oct 16, 2026

pelesh reviewed Oct 17, 2026

View reviewed changes

lukelowry added 5 commits October 18, 2026 02:42

Read OPF cases from JSON and rebuild components on ports

680fdde

Adopt the standard model interface and start from flat

8e91d10

Tag inequality constraints, store multipliers, and simplify parsing

f679ffd

Remove island checks

fe85267

distill documentaiton

5bc7404

lukelowry force-pushed the lukel/OPF-NEW-dev branch from 38dc4b5 to 5bc7404 Compare October 18, 2026 07:54

remove edge case

657ac9f

lukelowry requested a review from pelesh October 18, 2026 13:40

		/// Relative tolerance of the optimal cost against a reference cost
		constexpr double COST_TOLERANCE = 1.0e-8;

		* @tparam ModelT - model type with `VARIABLE_SIZE` local variables and
		* `CONSTRAINT_SIZE` local constraints

	[[maybe_unused]] __attribute__((always_inline)) static void mapped_store(double val,
	[[maybe_unused, gnu::always_inline]] static void mapped_store(double val,

		static constexpr RealT DERIVATIVE_TOL = 8 * std::numeric_limits<RealT>::epsilon();
		static constexpr RealT BALANCE_TOL = 1.0e-12;

		virtual int evaluateConstraints() = 0;
		virtual int evaluateJacobian() = 0;

Conversation

lukelowry commented Sep 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Checklist

Further comments

Uh oh!

pelesh commented Sep 28, 2026

Uh oh!

lukelowry commented Sep 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pelesh left a comment

Choose a reason for hiding this comment

Uh oh!

superwhiskers left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pelesh Oct 16, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nkoukpaizan left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pelesh left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lukelowry commented Sep 28, 2026 •

edited

Loading

lukelowry commented Sep 28, 2026 •

edited

Loading

pelesh Oct 16, 2026 •

edited

Loading